#include "../include/Distributions.hpp"
#include "../include/Density.hpp"
#include "../include/Misc/Laplace.hpp"
#include "../include/Kernel.hpp"
#include "../include/Particle.hpp"
//#include "../include/Resample.hpp"
#include "../include/utils.hpp"
#include "../include/Models/Likelihood.hpp"
#include "../include/Misc/EP.hpp"
#include "../include/Gibbs/TransdimGibbs.hpp"
#include "../include/Gibbs/TransdimGibbsst.hpp"
#include "../include/Gibbs.hpp"
#include "../include/Gibbsst.hpp"
#include "../include/Transdim_MCMC.hpp"
#include "../include/Misc/Poorman.hpp"
#define NU 200
#define _M 1000
extern "C"{
void Gibbs(int *M, int *n, int *p, double *Y, double *pX,double *res,double *res2)
{
	mat X=vect2mat<double>(pX,n[0],p[0]);
	mat bar(p[0],p[0]);
	bar.fill(10);
	mat S=diagmat(bar);
	//cin >> b;	
	Transdim_MCMC TdG(S,X,Y,(int)p[0]/2);
//	TransdimGibbs TdG(S,X,Y,(int)(p[0]/2));

	
//	Gibbss TdG(10,X,Y,1);
	mat G=TdG.Sample(M[0]);
	//vector<vector<int> > R=GreedySearch(G(span((int)(0.1*M[0]),M[0]-1),span::all),0.2);
	//vect2vect(R[0],res2);		
	mat2vect<double>(G,M[0],p[0]-1,res);

}
}
